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I. INTRODUCTION 


This thesis models numerically, the effects of a pressure disturbance incident on 
a fluid loaded reinforced plate using finite difference methods. The program produced 
explores the reaction of a plane wave striking the plate (hull of a ship) at different 
frequencies, for different plate lengths and for various spacings of the reinforcing beams 
(or members). 

Figure 1 shows the problem being modeled. The figure shows a cross section of 
a reinforced flat plate with water on the top side of it and a vacuum below it. For this 
problem, it is assumed that the plate is infinite in the x direction along the horizontal 
axis. Since the plate structure is periodic in space, the model will only look at one 
section of this infinite plate. The behavior at other sections is assumed to be the same. 
We will use three fundamental equations in the model. The first equation is the wave 
equation, the second is the linear-inviscid-force equation, and the third is the equation 
of motion for a thin plate. The coupling of these equations at the fluid-plate interface 
makes the problem very difficult to solve analytically. 

Chapter II derives the relationships among these three equations as they are used 
in the model. Assumptions used in the program are spelled out in this chapter. Chapter 
II also includes the steps necessary to scale each equation to facilitate its use in the 


computer program. 


Chapter III explains the operation of the FORTRAN computer program used to 
model the problem. The program takes as an input the incident and specularly-reflected 
sound waves and determines the scattered pressure due to the fluid-structure coupling. 
To make best use of the Naval Postgraduate School’s mainframe computer the main 


program is divided into a number of smaller subroutines. These are also described in 


Chaptér III. 


Chapter IV discusses the results of the program runs. It also reviews the 


weaknesses in the model and model limitations. 
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Figure 1. Cross Section of Hull 





Il. THEORY 


A. THE WAVE EQUATION 

It is known that sound travels through seawater, or any fluid medium, as a pressure 
wave. This is possible because seawater is a compressible fluid which is affected by 
pressure disturbances. To approximate the movement of sound in a fluid medium we 
have used the linearized, lossless wave equation for pressure as derived in Kinsler, Frey, 


Coppens and Sanders [Ref. 1:p. 104], 
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where ¢, is the fluid sound speed and p' is the total pressure. 

The pressure wave can be split up into various parts, each of which is a solution 
of the linear wave equation. In this case we want to be able to solve for pressures that 
would be present in the fluid which modify the pressure resulting from specular reflection 
of an incident plane wave from an acoustically hard surface. This known specularly 
reflected pressure wave is combined with the computed scattered pressure wave to make 
up the total reflected pressure wave. 


Symbolically, we have: 


p*=pl+p *+p*, C4 


where p’ is the total pressure, p' is the incident pressure, p® is the specularly reflected 
pressure, and p° is the scattered pressure. The total radiated pressure is a superposition 
of the reflected and scattered pressures. 
To eliminate the specularly reflected wave from the equations of motion we apply 
the boundary condition 
a! a) 2.3) 
dy ay 
In this case the normal derivatives of the incident and specularly reflected 
waves must cancel at the surface of the plate. The vertical coordinate Y is equal to zero 
at the surface. 
We assume that the source of the incident pressure is far from the surface so that 
a plane wave approximation is adequate to model the incident pressure, p’. 


The time harmonic incident and specularly reflected plane waves are of the form: 


x-d! , 
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where 
a cos0, (2.6) 
sin, 
and 
cos0, (2.7) 
sinO, 


are the unit direction vectors of the two plane waves. 

Both the incident wave, Equation 2.4 and the reflected wave, Equation 2.5 must 
satisfy Equation 2.1. Therefore, because of superposition, p* must also satisfy the wave 
equation, 
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B. THE LINEAR INVISCID FORCE EQUATION 
The linear inviscid force equation, taken from Kinsler, Frey, Coppens and Sanders 


[Ref. 1:p. 104], is basically the inviscid Navier Stokes equation in which nonlinear terms 


have been neglected. It is valid for acoustic disturbances of small amplitude, and can be 


written 


a Weer (2.9) 


where p, is the constant equilibrium density of the fluid, and w is the velocity of the 
fluid. 

Since we assume no cavitation of the fluid or penetration of the fluid at the fluid- 
plate interface, the component of the fluid velocity normal to the interface and the plate 
transverse velocity are equal at the fluid-plate interface at all times t. Letting w represent 
the transverse displacement of the plate (y direction), we have 


o : T 
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Through use of Equations 2.2 and 2.3 the total pressure p’ can be replaced by the 
scattered pressure p°: 
iw : S$ 
po = cae (2.11) 
of ey 
C. EQUATION OF MOTION FOR A THIN PLATE 
To properly understand what is happening on and within the thin plate, a reference 
system must be assigned. In Figure 2 we see that the vertical plane will be designated 


the y-axis in which the thin plate experiences transverse displacement w. The horizontal 


plane will be the x-axis which has longitudinal displacement u. To set up our equations, 


we assume that the centerline (axis) of the beam stays perpendicular to its cross section 


during bending. 





Figure 2 Plate Reference System 


As put forth by Kinsler, Frey, Coppens and Sanders [Ref. 1:p. 58], Hooke’s law 


states that if the strain is small the stress is proportional to it. The stress on the plate is 


described by the equation 


E &w 





where: 


1-v? dx? 


(Quid) 


0,, is the stress on the plate, E is the Young’s modulus, and vp is the Poisson’s ratio. 


Consider a small section of the plate of area AxAz with moments and forces applied 


as shown in Figure 3: 
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BLOW-UP OF SECTION A-AA 





Figure 3. Forces and Moments 
Note that there are no forces or moments applied which have z dependence. The 
motion therefore is assumed to be two-dimensional. Since the plate is not rotating, a 


force balance is performed on the moments and moment arm, neglecting terms of order 


Ax’, giving us 


OM _ (2.13) 
Ox 


where: 

M = bending moment on the element, 

q = load on the plate element (uniform for a normally incident wave), 
and F = flexural force applied to the element. 


Equating forces in the y-direction applied to the section of plate yields 


po ame Oe oA chr, (2.14) 
or Ox 


where: 

S = cross sectional area of the plate, 
p, = density of beam section, 

h = plate thickness. 


Dividing out the AxAz terms yields 





acai ae (2.15) 
Ox ot? 
Using the definition of the moment 
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where the moment of inertia is found from 


h° (2.17) 
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we obtain the equation for deflection of a plate 


BE ow ppcaliaoaet (2.18) 
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The loading of the thin plate is due to pressures exerted by a fluid in contact with 


the thin plate. Rewriting Equation 2.18 including the pressure terms yields: 





pa a (2.19) 
ax4 ot? 
in which 
3 
Eh (2.20) 
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is the flexural rigidity and h is the plate thickness, Junger and Feit [Ref. 2:p. 194]. 


D. THE SCALED EQUATIONS 


1. The Scaling Factors 


In order to use the three governing Equations 2.1, 2.9 and 2.19 in the finite 
difference code, they must be scaled or nondimensionalized. The steps are as follows. 


The independent variables x and t are scaled through use of the fundamental length L 


and the wave frequency w, respectively. 


x=—, and > EP 


ty | & 


where x and y are the scaled spatial coordinates, and L is the one half the distance 
between stiffeners measured in the x-direction. 


The equation for the scaled time is 


ap (2.22) 


where w is the frequency of excitation of the incident plane wave. 


We then have 
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Bc a (2.23) 
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and 
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Additionally, the dependent variables w and p must be scaled as follows. For 


the displacement of the plate (in the y-direction) we put 


w=, (2.25) 
W 


where w_ is the scaled vertical displacement of the plate and W is the displacement 


scaling factor. For pressure, put 


1] 
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where p is the scaled pressure and P, = Lup, W is the pressure scaling factor. 
The wave numbers for waves propagating in the fluid, and along the plate are 


given as follows, 





k=, (2.27) 
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respectively. 
Additional parameters resulting from the scaling analysis are 
ky 2.29 
joo a) 
aie 


where () is the squared ratio of fluid to plate wave numbers and @,, the coincidence 


frequency of the plate, is defined by 





p hey ; (2.30 
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We also define the dimensionless quantity 


peer (2.31) 
p JW, 


where ¢ is a dimensionless quality which is a measure of the fluid loading on the plate. 
2. The Wave Equation 


The wave equation governing the scattered fluid pressure was found to be: 


) he 
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Scaling the above equation yields: 
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which reduces to the desired non-dimensional form, 
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3. The Linear Inviscid Force Equation 


The linear inviscid force equation 


ew ap? 
Dye (2.36) 
at oy 
in scaled form reduces to 
Wer= De. (2.37) 


4. Thin Plate Equation 


The dimensional equation of flexural rigidity is: 
po +p Zep (2.38) 
Gun ot? 


By substituting in the scaling factors it becomes: 


Dvw----+p Aw L*Ww-- = -p "L?w*p w, (2.32) 


XXXX 


or, after more simplification, it becomes 


Wesszth, L’we: = Poe G95) (2.40) 

This is the scaled form of the thin plate equation that will be used in the 
computer program. Note that at y = 0, p® equals p'. Hence, p* has been replaced with 
p + 2p! 

5. Stiffeners 

In the present analysis the stiffeners will be treated in an idealized way. The 
end conditions we will assume are imposed at the point where the stiffeners are located 
will be equivalent to those of a clamped plate for which the boundary conditions are 


w = Oat x = +1, andw, =0 x = #1. 


E. THE FLUID DOMAIN 
As stated before, it is assumed the plate on which the sound wave impinges 1s 
infinite. That assumption makes it possible to establish the following conditions. The 


boundary conditions on the section of plate we are looking at will be periodic. That is 
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to say, for normal incidence of p', displacement on the plate at x = -1 will be the same 
as atx = 1. For arbitrary angles of incidence, the pressures atx = +1 differ only by 
a scaling factor determined by the angle of incidence. Figure 4 and the following 


equations will show how the time harmonic solution is dealt with in the problem. 
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Figure 4. Side View of Fluid Domain 





To understand how the fluid pressure is disturbed, the wave equation is considered 
for which a prescribed pressure is applied from the fluid to the surface of the plate. This 
uncoupled problem can be solved analytically by the method of separation of variables. 


The scaled wave equation is 


PtP aemne (2.41) 


ee tt 


Applying periodic boundary conditions results in 


p(-l,y)=p(,y) and p,(-1,y)=p,(1,y). (2.42) 


The forced pressure at the plate surface is given by the boundary conditions 


p(x,0,t)=g(x)e* = where g(-1)=g(1) and g(-1)=g,(1), (2-43) 
where: 
g(x) = prescribed pressure applied to the fluid by the plate. 
Assuming a time harmonic solution with equivalent frequency dependence to that 


of the forcing function, the scattered pressure can be expressed as 


PR Y.1)=o@ Ne C- 
By substitution of Equation 2.44 into Equation 2.41 we obtain the reduced 
Hehmholtz equation for (x,y), 
batdoscta°p =0, (2.45) 
with boundary conditions: 
o(-1y)=o0,y) and $(-Ly)=$,(Ly) (2.46) 


and 


(x, 0) = g(x). (2.47) 


Henceforth we will be dropping the bars from all equations, giving us the 


following, 
o,, 40,4 0° =0, (2.48) 
o(-1,y)=$(1,y), (2.49) 
o,(-1,y)=9$,(y), (2.50) 
(2.51) 


o,(-1,y) =o,C.y). 
To solve for ¢ the standard method of "Separation of Variables" is applied to the 


equation. First we write ¢ as a product of functions of one variable each, 


 =X(x)¥(y). (2.52) 


Upon substitution into the Helmholtz equation and separation of variables, we find a 
solution 
Nei (@.)= Cue ae (2.53) 


When solved for Y, two answers result: 


Y (y)=Aeove"'™ when a >nn, (2.54) 


and 


PAG) Ae “yvnint-a* when a<nNT. (2.555 


These solutions are chosen to avoid the physically implausible solution that the scattered 
pressure aS y approaches infinity will blow up represent an incoming wave. 
Combining our solutions we have 


= N oo 
o=-)> o(xy)=>> oe AEN aaeP eines ys One ane nena? (2.56) 
n=0 


n=0 n=N+] 
Note that as y approaches infinity the decaying terms become exponentially small. Now 


from our boundary condition at y = 0, a, can be found such that 
g(x)=)> ene with a,=5 fae “UARX Ty (2.399 
n=0 1 


It is important to note that in the numerical solution of the fully coupled problem, 
the radiation waves must be properly modeled as y approaches infinity. The numerical 
domain cannot be semi-infinite in extent in the y direction. It must be truncated at some 
value of y, which we call yg, where a boundary condition must be imposed to ensure 


proper modeling of the radiated pressure. 


F. THE PERIODICITY PROBLEM 

To ensure that a wave with any angle of incidence can be correctly modeled by the 
computer program the periodicity must be worked out. If we were to have a wave 
normally incident on the plate we would expect a solution for the pressure which is 


periodic in x to be 


b= 6,0)". a) 
n=0 


To account for an arbitrary angle of incidence, Equation 2.58 needs to be modified 
in the following manner. A correction factor which takes into account the phase change 
at x = -l vice x = 1 must be added. The modified equation for the scattered pressure 


1S. 


b=0F 6,0)e he Home, 2.59) 
n=0 


where 6, is the angle of propagation of the incident pressure wave. 


Continuing with this argument one can write 
7, 7~ kcosO,+nn. (2.60) 


Substitution of this x dependence into the Helmholtz equation yields 


d*o. 


dy? 
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Substituting these relationships back into our equation for ¢ gives us 


N : : pet oe ; 
b= ae PT eT ae, (2.62) 
n=-M —s0 N+} 


where 


Cae 
5 =) ee (2.63) 
nr ly? 42 k<y_, 
and in which the first sum includes all values of n for which k > y, and the latter two 


sums account for the remaining modes. 


G. FINITE-DIFFERENCE MODELING 
1. The Taylor Series 
Refer to Numerical Solution of Partial Differential Equations: Finite 
Difference Methods, Chapter 1, by G. D. Smith [Ref. 3] for a more complete explanation 
of finite difference methods. 
All three governing equations will be center differenced with truncation errors 
of order (Ax’, At’). 


Finite difference approximations for u,, U,,, Ux.x, are 


(u,),- 1 + 0(n), (2.64) 
U;,,~2U;+U,_, 3 ) 
(ie ae ) (2.65) 


and 


(U;,9 74,2) -4(u,,, +4,_,) +6(U)) 


74 +O(h 2 (2.65) 
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where u; = u(x). 
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2. Differencing of the Wave Equation 


We have already shown that the wave equation is of the following form 


Pee =e Alee., (2.67) 


xXx ay, 


and that we have the following conditions, 
x, = 1Ax, may t, = nAt, (2.68) 


where i, j, and n are integers. This allows us to describe the pressure as follows, 


P= P(x,y,t,) and  Ax=Ay=h. (2.69) 


So by substituting into Equation 2.66 we arrive at the final form, 


n+] 
a = 


ae 
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3. Finite Differencing of the Linear Inviscid Force Equation 
We have already shown that the linear inviscid force equation is of the 
following form 
ee, (2.71) 
oy 
This equation is applied at y = 0 or when j = O in the numerical domain. 
The displacement w does not depend upon y, so we introduce W" = W(x,, t, ) and write 


the differenced form of the equation as: 


mal 


wr? = =Ar7 i a som (2992) 
y) 


This equation is solved for the P®,_, which is exterior to the numerical domain. 
It can then be substituted into the differenced form of the wave equation to solve for 
p**!., , the pressure on the plate fluid interface at time level n+1 solving for P*,., we 


have 


p® = 289 yl oy "iy ae (2.73) 
| 42 { { { i,1 


A 
Note that W°*!; is first needed to solve for P",., which means a prescribed order must be 


followed at the fluid plate interface. First W°*'; is found, the P®;.,, and finally P"*?; o. 


4. Finite Differencing the Plate Equation 
The last of the major equations that we must put in the form of finite 


differences is the thin plate equation. The scaled thin plate equation is: 


XXXX 


~-Liekk* 
W—+k L‘W-- = Ee (ps 2P'), C 


Differencing this equation results in an expression for the displacement W at 
time level n+1 in terms of the preceding two time levels and the incident plane scattered 


pressure at time level n. 
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5. Finite Differencing Side Boundaries 
When differencing the wave equation at the side boundaries for the pressure 
we have a similar problem to that found in applying the wave equation at the fluid plate 
interface. The problem is that there are values of the pressure needed which are exterior 
to the numerical domain. Periodicity will be used as an additional condition to eliminate 


this problem. We have the scattered pressure in the form: 


PY b, (ene er: (2.76) 
n=0 


What we are doing here is taking the pressure at the node i = M+1, x = Ax 


and replacing the pressure 1 = -M+1 and x = -1+Ax. From Equation 2.76 we obtain 


Pas D dylyptemetrae Harare 2.79 
+1y oon 
Pragte emer) $0 p,)( = 1)"e int Ax, ikcos6, (2 . 78) 
and we also have: 
Ppa geSs iyi ventelene HV ane 2.79) 
-M+ nJ pin 


Now we can make the following two statements: 


ZS 


-ikA 8, -ikcos@ 
=e ikAxcos c ikcos DS b (-1remn* 


-ikA 8, ikcosé 
pense ikAxcos iol COs S o (-1)"ei"™4* 


Combining the previous statements we get: 


-2ikcos8 


Ro pn 
Pus j=e oy 


An equivalent argument used at the left boundary yields: 


n ___ 2ikcos8 
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6. Finite Differencing the Radiated Boundary. 


(2.80) 


(2.81) 


(2.82) 


(2583) 


The final edge of the numerical domain 1s the artificial boundary which must 


allow free passage of radiated pressure disturbances. It must be far enough away so that 


the evanescent modes can be neglected. To see how the boundary condition is derived 


we Start with the scattered pressure p* and try to derive an analytical expression for the 


pressures normal direction at the boundary, as is done in Scandrett and Kriegsmann 


[Ref. 4]. This type of radiation boundary condition 1s often referred to as a non-local 


radiation boundary condition. 


From previously we have 


N 


p*(xy,t) =e" { S- oem “Y"'+ evanesent modes), 


n=-M 


where 
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(2.84) 


Y,= ~ikcos®,+nn, (2.85) 
B = k2-y? (2.86) 


and 
Vp S=k2p 5. (2.87) 
Additionally in our problem we have 


k=kL and ke, (2.88) 


but as y approaches infinity p° is equivalent to only the radiated modes. We include the 


time dependence into the a, coefficient such that, 


N 
ps~ > ae, 2.89) 


n=-M 


and in the following steps use orthogonality to find a,(t)e”*"”. We have 


I nm . ] e e 
fe “Wet S(x,y,t)dx = S- a,(e"” fe tae rina (2.90) 
-] n=M =] 


and 


2 Vx", i OMY a Yal 
a, ix{ -ikcos8,+nn -(-ikcos8,+mr)] (Z 9 1) 
e ix[nn -mtr] 


We can therefore write 


where € = dummy variable of integration. 
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] 
ae =— fe =p S(Eynde, (2.92) 
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Noting the similarity to zero in the following equation 


N 





YS 2 +p, Shae elo, (2.93) 
n=-M 
Which can be rewritten as 
1 
op* 1 iy,(x-&) Op 5 
ea: ; Nd ae Q. (2.94) 
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Remember, these equations pertain to the artificial boundary which is modeling infinity. 
For this problem y is equal to y, at the artificial boundary. 
When we start to finite difference Equation 2.94, the artificial boundary and 


use the trapezoidal rule of integration we obtain: 


Pirye~Piny-1 te = iy,(%;-8) 7 pp mtl_ pm-ly_ _ (2.95) 
a ee (Puy ~Pry ))- =0 
2h 4Athe-M l=-1X 
in which we take: 
a l=+IX 7.96 
6, = 2 ( ‘ ) 
1 [#sIX. 


Rewriting Equation 2.95 we have 
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2h way? 81 Ber Dy(pret_prm}).0 (2.97) 
n=-M 
Or 
P ' De ; _ ; 
pions 3 Anan -Piiy )~ 0 (2.98) 
l=-IX 


Where the square matrix A is defined as 


A..= 0 Va 8) (2.99) 
an = 55 8, 


Combining the above equations gives us: 


IX 
Ee ee > A,, (Priy -Pi'y ) ~0 (2.100) 
1=-1X 
From the differenced form of the wave equation applied at j = IY (the artificial 


boundary) we have: 


Nt n nh 
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Pi = -P; a inn Pita Pine eae Oe 7242 )P., 100) 


Pity 


Combining Equations 2.100 and 2.101 and eliminating P"; ,,, we obtain: 


oy) 


nel _ nat ee n n n 
Piry = Pity Merrie cites ci) 


(2.102) 
at yA (Py Pr aioe 4 At* pr 
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Now we need to put this equation into a form which makes it possible for us 


to see the matrixes needed to solve the problem. First we write: 


pi alee 
[ae BA Puy 





mee 
1 Are -1 2.103 
Ph We 4 Pun * 93 [Pisvay* Pi-tayt2P iy] + (A 
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We define 
oan / te (2.104) 
2 
A= 8 ytd, (2.105) 
De ne 
and 
2 
B,=-3,t—— A), (2.106) 
te * 
Equation 2.102 can then be written in matrix form as 
APS Se eee (2.107) 
where 
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2 oe 
ee 2 4 Py (2.108) 


i h*k2 hk? 
Then the equation for the pressure vector at the next time level is given by: 
Pr’) = ABP" +C] (2.109) 
The pressure at P"*'. ,, is then assigned by taking the i* component of the 


vector ptt In solving the matrix Equation 2.107 we use the LINPACK library 


Gaussian elimination programs (CGESO and CGESL). 
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I. DESCRIPTION OF PROGRAM 


A. GENERAL 

The computer program written for this thesis simulates the effects of a pressure 
disturbance on a fluid loaded reinforced thin plate. A detailed description of the theory 
behind the program can be found in Chapter II. A pressure wave in a fluid medium is 
moving towards a reinforced plate which ts flat, and is supported by stiffeners which are 
at a constant spacing on the non-fluid side of the plate. The non-fluid side is assumed 
to be a vacuum. When the pressure sound wave hits the plate there will be some (a very 
small amount) distortion of the plate in the vertical direction. This displacement will 
Cause a disturbance in the fluid medium. The scalded pressure amplitude for different 
frequencies and beam spacing will then be plotted and the magnitude of the propagating 
modes will be recorded. 

The program stores the time dependent wave equation with time harmonic forcing 
by allowing transient effects to radiate away from the numerical domain. When a steady; 


state has ben achieved the code automatically stops the calculations. 


B. SUBROUTINES 
To make effective use of the Postgraduate School’s mainframe computer, the 
program has been broken down into a number of smaller programs (i.e. subroutines) 


which can be more efficiently programmed and run. 


30 


The main program is used to coordinate the subroutines and to pass information 


from one subroutine to another. The following are descriptions of the subroutines. 


1. Subroutine INIT 

The INIT subroutine initializes the computer program. In this subroutine we 
set up our two main matrices. The first of these is the "P" matrix, which represents the 
fluid pressure in the fluid domain. The second is the "W" matrix, this represents the 
displacement of the plate surface which bounds the fluid medium. 

Although we have the ability to change any number of parameters in the 
program, we will concentrate on two. These are the beam spacing, and the frequency 
of the pressure wave impinging on the plate. Both of these values are assigned in this 
subroutine. The INIT subroutine sets all of our counters and subscripts to zero at this 


time and initializes the matrices which will be needed to model the artificial boundary. 


2. Subroutine DOMAIN 
The DOMAIN subroutine calculates the pressure at all nodes in the interior 
of the fluid medium at a subsequent time level in terms of the preceding pressure values. 
The subroutine calculates the pressure at the next time step by use of the differenced 
wave Equation 2.70. 
The equation above must be modified to model what is happening at the side 
boundaries of the fluid domain. As was stated earlier, our model is looking at a small 
section of an "infinite" fluid plate interface. This allows us to use the pressure in the 


next section to model what is happening in the section in which we are interested. This 


el 


1s necessary since otherwise we would have to use "undefined" points in our differenced 


wave equation at the edges of the fluid domain. 


3. Subroutine FLPL 

The FLPL subroutine models the fluid-plate interface along the top of the 
reinforced thin plate. Again, with finite differencing and the assumption that we are only 
looking at a small section of an infinite interface, we are able to model the effects of the 
pressure wave along the entire fluid-plate interface. This subroutine is also where we 
introduce our driving function into the computer program. This function is what starts 
things moving towards a steady state. Although this is not the most complicated 
subroutine in the entire program it 1s the one which has proved to be the hardest to make 


work as predicted. 


4. Subroutine ARTF 
The ARTF subroutine is the most complicated of any subroutine in the 
computer program. It models the artificial boundary of the fluid medium at infinity. 
Again, this is accomplished with finite differencing and matrix operations. 
5. Subroutine ALPHA 
The ALPHA subroutine calculates the magnitude of the radiating modes of the 
pressure wave. This data is gathered in the results section of this study. 
6. Subroutines SHIF and FINN 
The SHIF subroutine moves the entire problem to the next time level. This 


is done after the other subroutines have completed all the required calculations needed 
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to satisfy our domain and all of our boundary serine After exiting the SHIF 
subroutine the entire program is run again with the exception of the INIT subroutine. 

The FINN subroutine is a check to make sure that we have reached a steady state 
situation before we take any measurements of the plate displacement or of any fluid 
disturbances. Figure 5 is a graphical representation of the displacement convergence 
factor approaching zero. This was arrived at with a assigned frequency of 3500 Hz and 


L = 1.0. 


7. Subroutine CLOSE 
The CLOSE subroutine is used to send the data generated in other parts of the 


main program to the plotting devices. 


C. PROGRAM VERIFICATION 

Verification of the program was accomplished by the following method. We solved 
the waveguide problem (section 2E) for a known value of g(x). We let g(x) = sin ax and 
e'™*. This was accomplished by replacing the FLPT subroutine with a subroutine named 
TEST. As before the program was run until a steady state solution was reached. These 
results were checked against analytical solutions and found to compare very well. For 
g(x) = sin wx we graphically reproduced the expected results. For g(x) = e'™* we 
calculated the coefficients a, and found a, ~ | and all other a,’s ~ 0. The small error 
found in the calculation can be attributed to truncation errors in the finite difference 


approximations. These results match with analytical calculations. 
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Figure 5. Convergence Graph 


34 





IV RESULTS AND CONCLUSIONS 


A. INITIAL CONDITIONS 


To properly apply the computer program generated, inputted data, or conditions, 
must be within the tolerances of the main equations used. This 1s demonstrated by the 
use of .01 meters as the thickness of the thin plate and a frequency range from 500 Hz 
to 4000 Hz. Other inputted data, including the densities of steel and water were taken 
directly from Kinsler, Frey, Coppens and Sanders [Ref. 1:pp. 461-463]. 

The choice of the Ax and Ay terms is based on the need to have a significant 
number of data points in the waveguide domain that will give an accurate picture of what 
is happening. This is also true for the time step chosen, we need to make sure that the 
At was small enough to keep the problem from blowing up. A Courant-Lewy-Friedrichs 
Stability condition must be satisfied for the domain which limits the size of At. 
Additionally as was mentioned before, the computer program was run for the same 
amount of time for each trial, this time period is long enough to ensure that a steady state 


solution was being achieved. 
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B. PROGRAM OUTPUT DATA 

The data obtained from the computer program demonstrates the following major 
points. When the frequency is increased you get more propagating modes from each 
particular plate length. When L is changed by a specific amount you do see a difference 
in the magnitude and shape of the deformation graphs. 

Our first table shows the calculated magnitudes of the propagating modes of the 
pressure wave at an assigned frequency of 750 Hz at different beam spacings. 


Table I. MAGNITUDE OF RADIATING MODES. 


For a frequency of 750 Hz 
with L = 0.8 ay = .467808127 


with L = 1.0 ay = .484098911 
a, = .307484090 


with L = 1.2 a» = .232444227 
a, = .367214561 


In addition to the data in Table I, we have in Figures 6 through 8 a graphical 
representation of the scalded pressure amplitude. The pressure was measured at the 
artificial boundary , y, at a frequency of 750 Hz. The horizontal axis is of a length of 
2L for each of the graphs starting at x = -1 and going to x = 1. The numbering on the 


horizontal axis represents each of the 41 pressure nodes used in the finite difference 


code. 
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SCALED PRESSURE AMPLITUDE 


L=0.8 FREQUENCY = 750 HZ 
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Figure 6. Sample Pressure Graph 
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SCALED PRESSURE AMPLITUDE 


L= 1.0 FREQUENCY = 750 HZ 
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Figure 7. Sample Pressure Graph 
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SCALED PRESSURE AMPLITUDE 


L = 1.2 FREQUENCY = 750 HZ 
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Figure 8. Sample Pressure Graph 
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Table II is a list of the magnitudes of the a,’s for program runs with a frequency 
of 3500 Hz. 


Table 1. MAGNITUDE OF RADIATING MODES. 


For a frequency of 3500 Hz 


with L = 0.8 = ap = .451181054 
a, = .323620617 
a, = .088650584 
with L= 1.0 = ag = .365310729 
a, = .260265648 
a, = .062102031 
a, = .060141488 
a, = .076183497 


with L = 1.2. a = .252027571 
a, = .165709317 
a, = .047427550 
a, = .046469088 
a, = 049672558 
a, = .043821129 
a, = .067380726 


From the data above we observe that at a frequency of 3500 Hz there is a 
substantial change in the magnitude of the propagating modes as L 1s increased from 0.8 
to 1.2. It can be seen from other runs that the higher the frequency the larger the 
differences in the magnitudes of the propagating modes. 

Figures 8 through 10 show the scaled pressure amplitude data obtained from the 
program runs with an assigned frequency of 3500 Hz. It is observed that changing of 
the length L on which the Ax is scaled has an effect on not only the amplitude of the 


pressure wave but also their shape. 
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SCALED PRESSURE AMPLITUDE 


L=0.8 FREQUENCY = 3500 HZ 
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Figure 9. Sample Pressure Graph 
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SCALED PRESSURE AMPLITUDE 


L = 1.0 FREQUENCY = 3500 HZ 
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Figure 10. Sample Pressure Graph 
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SCALED PRESSURE AMPLITUDE 


L=1.2 FREQUENCY = 3500 HZ 
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Figure 11. Sample Pressure Graph 
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Table III is a list smaple energy calculations taken at three frequencies inputted into 
the computer program. The data presented demonstrates that for theses frequencies and 
lengths run in the computer program that the fundamental mode contains the largest 
amount of energy. Also evident is that the energy contained in the fundamental mode 
starts becomes more constant as the frequency is increased for different lengths of L. 


This data is in agreement with the theory of sound propagation. 
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Table If ENERGY CALCULATIONS 





Percent of Energy in Percent of Energy in 
Propagating Mode Propagating Mode 
For a frequency of 750 Hz For a frequency of 3500 Hz 
with L = 0.8 E, = 100% with L=0.8 E, = 66% 
I ee ve 
payne — 10 E, = 99% E, = 2% 
E, = 1% 
mths 152 E, = 77% Withee — 1.0 E, = 65% 
E, = 2% 
For a frequency of 1750 Hz EF, = 1% 
E,= * 


with L = 0.8 E, = 77% 


E, = 23% with L = 1.2 E, = 64% 
E, = 27% 
with L = 1.0 E, = 68% E, = 2% 
E, = 31% E, = 2% 
E,= 1% E, = 2% 
E;= 1% 
with L = 1.2 E, = 64% E, = 2% 
E, = 32% 
E,= 3% 
E, = 1% 


* denotes much less then 1% 
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Table IV shows the magnitude of the primary propagating mode for each of the 
different frequencies and lengths inputted into the computer program. 


Table IV MAGNITUDES OF PRIMARY MODES 


Magnitude of a,’s 





Eas Dai) C7 

Frequency 

750 Hz 467808127 =—.484098911 = .367214561 
1000 Hz .6182562721 .444052100 = .262355387 
1250 Hz .553467214 457091272  .343211770 
1500 Hz .498957574 .415586054 .320898831 
!/S00Hz .404151559 431289613  .338643909 
2000 Hz .486700475 .433078885 = .3 16784620 
2250, 394520938 .400619507 = .252146780 
Za 0 oer .385663450 = 392564356 ~—.267684340 
2750 Hz .385663450 = .392564354 = .267684340 
3000 Hz .420025572 353338059 ae 
3250 Hz .369634628 .226448834  .194508076 
3500 Hz .451181054 365310729" SezZs 20277 
S750 tae .446322470 = .330691218 = .226533830 
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The data obtained from the computer program and represented in Table III shows 
that there is always a difference in the magnitude of the fundamental propagating mode 
for changes in L at all frequencies comprised within the computer runs. The difference 
in the magnitude of the primary propagating mode between the runs with a length "L" 
of 1.2 and the runs with a length L of | are greater then the difference between the runs 
with lengths of 0.8 and 1.0. In the frequency range covered when L is increased from 
1.0 to 1.2 there is always a drop in the magnitude of the fundamental propagating mode. 
This is not always true when L is decreased to 0.8 from 1.0. Although the magnitude 
of these increases are smaller then the magnitude of the decreases this demonstrates that 
by increasing the spacing of the reinforcing beams will not always result in a decrease 
in the amplitude of the fundamental propagating mode in this frequency range. From this 
result we can conclude that additional factors play an important role in determining the 
magnitude of the propagating modes. 

Because the model is allowed to run until a steady state solution is reached there 
iS no appreciable difference in the absolute pressure measured at the artificial boundary 
for any of the frequencies observed. However, it is evident that as the length of the plate 
"L" is increased there is a slight drop in the pressure measurement. This is true for all 
frequencies. Figure 12 isa sample plate deformation graph for L = 1.0 and a frequency 
of 1000 Hz. 

In conclusion, the computer program produces results as expected for the different 
frequencies of the plane wave and lengths of the thin plate. The superpositioning of the 


reflected plane wave has showed that as L is increased there is a larger number of 


47 


propagating modes accompanied by a change in the amplitude of each mode. Additional 
results could be explored for different angles of deflection as well as different plate 
thicknesses. This will require modification of the computer program as it stands to 


gather data for additional angles and thicknesses. 
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SCALED DISPLACEMENT OF PLATE 


SCALED LENGTH OF 1.0 
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Figure 12. Deformation Graph 
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APPENDIX: COMPUTER PROGRAM 


FILE: 280CT FORTRAN Al 


—_— 


CHE HE HE BE HE IE HE IE HE IE HE HE HE IE IE HE IE HE HE IE FE HE HE IE HE IE HE FE HE HE FE IE HE HEHE HE HE HE K HE HEHE HE HE HE HE HE HE HE HE HE HE HE HE HEH HEHE 


AFTER A SOUND WAVE HAS HIT 


PRESSURE CEXCESS) D 
HORIZONTIAL POSITION Y 

W DEFORMATION KF 
SPD WAVE SPEED KP 
FREQ FREQUENCY OMEGA 


COUNTERS: I J N 


THE VALUES IN THE PARAMETER 
PARAMETER CIX = 20, IY = 80) 
PARAMETER CH = .025, DELTAT = 
PARAMETER CNDIM = 2xIX+1) 


p 
X 


AANQIANAIANIAININNINOND 


» THIS PROGRAM SIMULATES THE DISTRUBANCE IN A FLUID 
A PLATE. 


FLEXURAL RIGIDITY 
VERTICAL POSITION 
FLUID WAVE NUMBER 
PLATE WAVE NUMBER 
RATIO OF CKF/KP ) 2 


STATEMENTS ARE CONSTANTS 


20525 > 


S THE FOLLOWING ARE COMPLEX ARRAYS 


COMPLEX*16 
COMPLEX 
COMPLEX 
COMPLEX 
COMPLEX 
COMPLEX 


PC -IX: IX, 
LB, FACR 

AC -1a 
AAC NDIM, 
BBC NDIM, NDIM) 
CC -IX: IX) 
COMPLEX ZC NDIM ) 
COMPLEX®16 FL 
COMMON/’SS/P, W 
COMMON/TT/KP,KF,INFA, 
l EISP, OMEGA, 
COMMON/CC/LB, FACR 
COMMON EE/WGTWO, WGTHE, 
COMMON/PT/PGTWO, PGTHE, 
COMMON/HH/MINN, MAXX 
COMMON/IIZAA, BB 
COMMON/KK/ZIPVT 
COMMON/MM/NN 

COMMON 00/C 
COMMON/FFZA 
COMMONZUPZAL 
COMMON/RN/RUNS 
COMMON/FEL/U,V,TP 
COMMON/SN/FREQ, SPD 
COMMON/ GF/WGGRF , PGGRF 


INTEGER RUNS, IPVTCNDIM), 
PARAMETER CSAM = 1000) 
REAL KF, KP, INFA, 
REAL ¥8 
REAL ¥8 


0: 


me 
NDIM) 


WGTOT, 
PGTOT 


EISP, 


PGTHO,PGTHE,PGTOT 
REAL WGGRFC SAM ), 
REAL THETAI, FREQ, 
REAL X, T 

PI = 3.141592654 

CON = .015D0 


SPD 


a 


IY, 
IX) 


THETAI,L, 


SAM 


OMEGA, 
WGTWO,WNGTHE,WNGTOT,CON 


3), WC -IX: ex 33 


X, PI, 


CON 


L,PI 


PGGRFC SAM ) 
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elelelexe 
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AO 
— © 


80 


767 
768 
196 
487 
735 
629 
919 


WGTWO = 0.0D0 
WGTHE = 0.0D0 
PGTWO = 0.0D0 
PGTHE = 0.0D0 


CALL XUFLOW (0) 
CALL INIT 
RUNS = 0 
CALL FLPL 
CALL DOMAIN 
RUNS = RUNS + 1 


CALL ARTF 
IF CRUNS .GT. 1999) THER 
GOTO 60 
ELSE 
GOTO 61 
ENDIF 
CALL FINN 
T = T + DELTAT 
CALL SHIF 
IF CRUNS .EQ. 3000) THEN 
GOTO 80 
ELSE 


CALL ALPHA 
CALL CLOSE 


WRITEC2, 767 )THETAI, KE 

FORMATC1X, 'THETAI = ',F12.6,3X,' KF = ',F12.6) 
WRITEC2, 768) OMEGA, EISP 

FORMATCIX, * OMEGA = °,F12.6,3X,'EISP = ',F1l2.7) 
WRITEC2,196)DELTAT, KP 

FORMATC1X, 'DELTAT = ',F12.6,3X," KP = ',F12.6) 
WRITEC2,4987)H,L 


FORMAT(C1X,° = *,F12.6,3X,° L = ',fF6.3) 
WRITEC2, 7 55) FREQ, SPD 

FORMATC1X, "FREQ = ',F12.6,3X,'SOUND SPEED = ',F12.6) 
WRITEC2,629)PI 

FORMATC1X," PI = ',F12.7) 

WRITEC2,919) 


FORMATC1X, "FACR 
WRITEC2,%) FACR 


ue) 


ap 


AQAAQ 


WRITEC2,222) 


Zee FORMATC1X, 'FINAL DN ea X mi 
8 ¢. 8 ) 


XX = IH 


W , 
DO 333 I = -IX, IX, 1 


353 WRITEC2,666)XX,CDABSCWCI,2)3),CDABSCPCI,IY,3)) 


666 FORMATC1X,5F10.5) 
WRITEC2,223)RUNS 

223° FORMATC1X,*RUNS = 
WRITEC2,678) 


WRITEC2?, XJWGTOT 
WRITEC2,*®PGTOT 
STOP 
END 


SUBROUTINE INIT 


OOOO ee 


PARAMETER CIX = 20, 
PARAMETER CH = .025, 


*,1I7) 
678 FORMATC1X, BELOW IS THE FINAL WGTOT, AND PGTOT*) 


DEL 


THIS SETS UP THE INITIAL CONDITIONS OF THE 
PROBLEM INCLUDING THE DOMAIN AND THE 
ARTIFICIAL BOUNDARY 


THE VALUES IN eee 


ees STATEMENTS ARE CONSTANTS 
0 
TAT = .0325) 


PARAMETER CNDIM = 2xIX+1) 
c THE FOLLOWING ARE COMPLEX ARRAYS 


COMPLEX*16 EC ~LxX 


IX, 


COMPLEX LB, FACR 
COMPLEX AC -IX:IX, 
COMPLEX AAC NDIM 
COMPLEX ZC NDIM 3) 


COMMON/SS/P, W 


COMMON/TT/KP,KF,INFA, 


FILE: 280CT FORTRAN Al 


O: ITY, 33, WC -IX: IX, 3) 


-IX: IX) 
» NDIM )}, BBC NDIM » NDIM 


THETAI,L, X, PI, 


] EISP, OMEGA, T 


COMMON’CC/LB, FACR 
COMMON/FF/A 


COMMON/’GG/GAMMA, BETTA 


COMMON/HH/MINN, MAXX 
COMMON/ZII/ZAA, BB 


INTEGER I, J, N, LJ, 


KF, EISP, OMEGA, TIN, L, KP 


XIL, RCOND 


COMMON/KK/IPVT 
COMMON/FEL/U,V,TP 

: COMMON’SN/FREQ, SPD 
REAL THETAI, INFA, 
REAL ART, K, XI, 
REAL GAMMAC-G0: 40), 
REAL SPD, FREQ 
REAL X, T 


NN, 


BETTA (-40:40), PI 


U, V 


TP 
INTEGER NMODMIN, NMODMAX, MINN, MAXX 


INTEGER IPVTCNDIM) 
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) 


QOaqq oO 





= Gi 

LB = CMPLX(0.0, 1.0) 
FACR = CMPLX(1.0, 0.9) 
SPD = 1500.0D0 

FREQ = 3750.0D0 

N = 


= 90 

I = 0 

J = 0 

NN = J 

[ sx 

THETAI = PI72.0D0 

T = DELTAT 

NMODMIN = 0 
NMODMAX = 0 
KF = (2D0XPIXFREQXL)/SPD 

. GT = €€195000000000 .0XC .05%*3))71201-C. 282) ) 

KP = CCCCFREQ¥Z¥PI)DX¥2)%7700XH)I/FR)I XX. 25 
KP = 40.7534975147D0 


OMEGA = CKF/KP) 2 
INFA = CDELTAT#X2)/(C CKF 2) CL XX2) XCHRR2 DD 
DO 12 I = -IX,IX,1 

DO 11 J = 0,1Y,1 


PCI,J,1) = CMPLX(0.0, 0.9) 
PCI,J,2) = CMPLX(0.0, 0.0) 
PCI,J,3) = CMPLX(0.0, 0.0) 
11 CONTINUE 
12 CONTINUE 
DO 21 I = -IX, IX, 1 
WCI,1) = CMPLX(0.0, 0.0) 
WOI,2) = CMPLX(0.0, 0.0) 
WOI,3) = 0.0) 


CMPLX(0.0, 
21 CONTINUE : 


THIS SETS UP THE ARTIFICAL BOUNDRAY 


K = KFXL 
TIN = -KXCOSCTHETATI) 
DO 15 N = -40, 40, 1 
GAMMACN) = TIN + PIN 
IF (© ABSCGAMMACN)) .LE. K) THEN 
NMODMIN = MINCN, NMODMIND) 


NMODMAX MAXCH, NMODMAX) 
ENDIF 
15 CONTINUE 
MINN = NMODMIN 
MAXX = NMODMAX 


DO 18 N = MINN, MAXX, 1 
BETTACN) = SQRTCCK®®2) - CGAMMACN)®X2)) 
18 CONTINUE 


DO 78 I = =e IX, i 
We 79 LJ = “1X, Ix, 1 
ACI,LJ) = CMPLX(C0.0,0.0) 
79 CONTINUE 
78 CONTINUE 


ART = (2®CH®¥2))/7(0GXDELTAT) 
DO 300 Te = 1X, 51X, 1 
DO 301 LJ = =X; IX, } 
XI = IH 
XIL = LJ*¥H 
IF CABSC(LJ) .EQ. IX) THEN 


53 


C 


MQAMAQ O 


DEL = .5 
ELSE 
DEL = 1 
ENDIF 
DO 302 N = MINN, MAXX, 1 
ACI,LJ) = ACI,LJ) + ARTXDEL 
X¥BETTACN) XCEXPCL BXGAMMACN)®C(XI-XIL)) 


302 CONTINUE 


301 
300 


421 
420 


CONTINUE 
CONTINUE 
DO 420 U = 1, NDIM, 1 
DO 421 V = 1, NDIN, 1 
T =°=1XtUe 
LJ = -IX+t+V-1 
AACU,V) = INFAXCACI,LJ)) 
BBCU,V) = INFAXCACI,LJ)) 
IF (I .E@. LJ) THEN 
AACU,V) = AACU,V) + 1.0 
BBCU,V) = BBCU,V) - 1.0 
ENDIF 
CONTINUE 
CONTINUE 


CALL CGECO CAA,NDIM,NDIM, IPVT, RCOND,2Z) 
RETURN 
END 


SUBROUTINE FLPL 


THIS SUBROUTINE DEFINES THE FLUID-PLATE INTERFACE 
THIS IS WHERE Y=0 ALONG THE TOP OF THE PLATE. 


THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER CIX = 20, IY = 80) 
PARAMETER CH = .025, DELTAT = .0325) 
PARAMETER CNDIM = 2x*IX+1) 
THE FOLLOWING ARE COMPLEX ARRAYS 
COMPLEXX16 PC -IX: IX, O: IY, 3), WC -IX: IX, 3) 
COMPLEX LB, FACR, PTEMP 
CONPLEXX16 EL 
COMMON’SS/P, W 


sepa | | ORC THETAI,L, X, PI, 


EISP, ONEGA, T 
COMMCH’”CC/LB, FACR 
COMMON/RN/RUNS 


REAL KP, L, OMEGA, EISP, KF, INFA, PI, THETAI 
REAL X, 


7 
REALX4 SQ, Sid, BW1, BW2, BH3, BMS 
INTEGER I, J, RUNS 


EISP = 0.13101864566D0 


J = 0 

SQ = CDELTAT/CHXKF) ) x2 

SH = CDELTAT/ CHXKP )XX2)xx2 
Bil = 2-6*SW 

BW2 = G*SW 

BW3 = EISPXDELTAT*X2*KF/OMEGA 
BWG = -SW 

A = FLOATCIX)/FLOATC2*IX) 
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QO 


oleke) 


eoleke) 


- 


1 
10 CONTINUE 


pei. = INFAXCPCI+1,J5,2)+PC 


eel ,J,3) = INFAXCPCI4+1,J,2)+FACREP(-I-1 
( 


ne? ,J,3) = INFAX ( FACRXP (- Ett, ZIP CI), J, 2)-SXP 


a 


THIS LOOP COMPUTES THE INTERIOR NODES 


DO 10 I = -IX + 2, IX -2, 1 
X = IH 


WCI,3) = BWI¥WCI,2) -WCI,1) + BWN2xXCWCI41,2)4tWCI-1,2)) 
1 


+BWG*XCWCI+2,2)04+WCI-2,2)) 

+BW3XCPCI,0,2)+ FICX,T)) 
FEM = -(2*XH/DELTATXX2)*CUCI 
(I,J,3) = INFAXCPCI+1,J,2)4+PCI-1 
PTEMP d+2¥PCI,J,2) 


,3)-2EWCI,2)4*WCI,1)) + 
pole aoe C1, 2)+PCI, J 
Sela sk) 


THESE LOOPS COMPUTE THE NODES AT -19 AND 19 


= ~IX+1 

= (-IX+1)*H 

C(I,3) = BWI¥WCI,2)-NCI,1) +BW2*¥CWCI+1,2)4+WCI-1,2)) 
+BWNGXCWCI+2,2)4+NCI, 23) 
+BW3¥CPCI,0,2)+ FICX,T)) 

PTEMP = —-(2*H/DELTAT#¥¥2)*CNCI 

PCs), 3) = INFAXCPCI+1,J,2)+PCI-1] 

FTENP d)+2ePCI,J,2) 


-2*WCI, 
2)-G¥*PC 
(I,J,1) 


, 


,3) 2)4+WCI,1)) 

rd, Iencore cl, 

-P 

= IX -1 

= (IX-1)*H 

WCI,3) = BHIXNCI,2) - WCI,1) + BW2*CWCI+1,2)+WCI-1,2)) 
+BHGXCWCI, 2)+NCI-2,2)) 
+BW3X(P(I,0,2)+ FICX,T) 

PTEMP = “CoRHCDELTAT #x2 )x( WC x1 


4 
J 


m NK 
= vw 


) 

I,3)-2¥"WCI, 

1,J,2)-4XPC 

EePeEMP d+2RPCI,J,2) -PCI,J,1) 
THESE ARE THE END POINTS 

eh akee ee COS STHETAT)) 


J = 0 
WCI, 3) 
PTEMP 


, 


CMPLX(0.0,0.0) 
P(I,1,2) 


pln ce 
I,J, 


or 


+PCI, J+1,2)+PTEMP )+2xP 
FACR=CEXP(- 2*L BXLXKFXCOS(THETAI)) 
i= 10 
WCI,3) 
PTEMP 


) 


CMPLX(0.0,0.0) 
P(I,1,2) 


(I,J, 
+PCI,J+1,2)+PTEMP )+2"P CI, i 2) =P Clas; 
RETURN 
END 
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P(I,1,2) 
+#1,2)+ 


+P(I,1,2) 
+1, 2)+ 


PLEE: 


Zt 


280CT FORTRAN Al 


SUBROUTINE DOMAIN 


THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER (IX = 20, IY = 80) 
PARAMETER CH = .025, DELTAT = .0325) 
PARAMETER (NDIM = 2X®IX+1) 

THE FOLLOWING ARE COMPLEX ARRAYS 
COMPLEX*16 PC -IX: IX, O: IY, 3), WO -IX: IX, 3) 
COMPLEX LB, FACR 
COMMON’SS/P, W 


ON ee: KF, INFA, THETAI, L, X, PI, 


EISP, OMEGA ,T 
COMMON/CC/LB, FACR 


REAL KP, KF, INFA, THETAI, L, PI, EISP, OMEGA 
REAL X, T 
INTEGER I,J 


DO 22 I = 1-IX, IX-1,1 
DO 21 J =il, IY-l, 
PCI,J,3) = INFAXCP 


ve = 


CONTINUE 


22 CONTINUE 


2 
41 CONTINUE 


2 
31 CONTINUE 


Ga 


jee = INFAXCPCI4+1,J,2)+ FACRXPC(-I-1 


LB = CMPLX (0.0, 1,0) 
LA = L¥KFXCOSCTHETAI) 
wee I CEXP((2)*CLBXLA) ) 
= -IX 
DO 41 J = l, iyo ? l 
»J,2) 
-G¥PCI,J,2)4+PC1,J3+1,2)4+P01,J3-1,2))+ 
2xECL; J, c)-F Cl 


LB = CMPLX (0.0, 1.0) 

RA = L¥KFXCOSC(THETAI) 

ee CEXP((-2) XC LBXRA) ) 

DO ot J = l, Lys), s 

PCI,J,3) = INFAXCFACREXPC(-I+1,J,2)+ PCI-Il, 
-~G*¥PCI,J,2)4+P CR 4+1,2)+P Cie 
ext Clee )-F Cl, le 


J,2) 
,2))+ 
RETURN 

END 
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CG 


QOQN0QM 


SUBROUTINE ARTF 


THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER (IX = 20, IY = 80) 
»PARAMETER CH = .025, DELTAT = .0325) 
PARAMETER (NDIM = 2%IX+1) 

THE FOLLOWING ARE COMPLEX ARRAYS 
COMPLEX*16 PC -IX: IX, 0: IY, 3), WC -IXk: IX, 3) 
COMPLEX LB, FACR 


COMPLEX AC -IX: IX, -1IX: 1X) 

COMPLEX AAC NDIM » NDIM ), BBC NDIM » NDIM 
COMPLEX BCNDIM) 

COMPLEX CC -IX: IX) 

COMPLEX ZC NDIM ) 


COMMON/SS/P, W 

COMMON/TT/KP,KF,INFA, THETAI,L,X, PI, 
EISP, OMEGA ,T 
COMMON/CC/LB, FACR 

COMMON/FF/ZA 

COMMON/II/ZAA, BB 


SUBROUTINE SHIF 
THIS MOVES THE PROGRAM TO THE NEXT TIME LEVEL 


THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER CIX = 20, IY = 80) 
PARAMETER CH = .025, DELTAT = .0325) 
PARAMETER CHDIM = 2x*IX+1) 
THE FOLLOWING ARE COMPLEX ARRAYS 
COMPLEX*16 PC -IX: IX, QO: IY, 3), WC -IX: IX, 3) 
COMPLEX LB, FACR 
COMMON/SS/P, W 
COMMON/TT/KP,KF,INFA, THETAI, L, X, PI, 
] EISP, OMEGA, 
COMMNON/CC/LB, FACR 


REAL KP, L, OMEGA, EISP, KF, INFA, PI, THETAI 
REAL X, T 


INTEGER I, J 


a 


) 


1 
G2 CONTINUE 


10 


30 
20 


150 


46 


ns 


COMMON/’FEL/U,V,TP 


COMMON/KK/IPVT 


COMMON/O00/C 

REAL KP, L, OMEGA, EISP, KF, INFA, PI, THETAI 
X, T 

ees IPVTCNDIM), JOB ,U,V, TP 

CMPLX (0.0, 1.90) 

L¥KFXCOSCTHETAT) 

L¥KFXCOSCTHETAT) 

DO 42 I = -IXt+l1l, IX-l, ] 

) 


CCI) = INFAXCPCIF+1, ITY,2)+PCI-1,1Y,2)+2xPCI,1Y-1,2)) 
+ C2 5 CGXINFAJ)XFPCI,IY,2) 


ir 
> 
aout ot 


I = -Ix 
FACR = CEXPC C2) *CLBXLA)) 


| “aa = INFAXCPCI+1,1Y,2)+FACREPC-I-1,1Y,2)+2PCI,1Y-1,2)) 


+ (2 - CGXINFA))XPCI,IY,2) 


I = Ix 

FACR CEXPC(-2) CL BXRA) ) 

CCT) INFAXCFACREXPCI-1,1Y,2)+PC-I+1,1Y,2)+2*PCI,1Y-1,2)) 
+ (2 - CGXINFADDXPCI,IY,2) 


DO 10 I = 1, NDIM, 1 
BCI) = INFAXCCMPLX(0.0, 0.0)) 
CONTINUE 


DO 20 I =1 
DO 30 LJ 
BCI) = B 

CONTINUE 

CONTINUE 


DO 150 [I = 1, NDIM, 1 

BCI) = BCI) + CC-IX+I-1]) 
CONTINUE 
CALL CGESLCAA,NDIM,NDIM, IPVT,B, JOB) 
DO 46 I = -IX, IX, 1 

PCI,IY,3) = BCIX+1+I) 
CONTINUE 
RETURN 
END 


pL JIKPC-IXFES-2 
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22 
19 


23 
20 


29 


40 


’ 


DO 19 I = -IX, IX, 1 
pO 22 J = 0, IY, 1 
P(t, v,1) = PCI,J,2) 
CONTINUE 

CONTINUE 


DO 20 I = -IX, IX, 1] 
DO 23 J = 0, IY, 1 
PCI,J,2) = PCI,J,3) 
CONTINUE 

CONTINUE 

DO 29 1 = -!I I 
WCI,1) = , 

CONTINUE 

DO 40 I = - 
WCI,2) = 

CONTINUE 

RETURN 

END 


FUNCTION F1CXP, TIME) 

PARAMETER CIX = 20, IY = 80) 

PARAMETER CH = .025, DELTAT = .0325) 

COMPLEX*®16 PC -IX: IX, O: IV, 3), WC -IX: IX, 3) 
COMPLEX*®16 Fl, Sl 

COMMON/SS/P, W 

COMMON/RNZRUNS 


I 


=i 


X, IX, 
WCT,2) 
xX x 
W > 


(1,3) 


COMMON/TT/KP,KF,INFA, THETAI,L, X, PI, 
] 


EISP, OMEGA, T 


REAL KP, L, OMEGA, EISP, KF, INFA, PI, THETAI 
REAL X, T 
REALX8 Q9,Q5,5 


IFCRUNS .EQ. 1) THEN 
S = 0.0. 
EESE 

S = TIME+XPXKFXCOSCTHETAI) 
ENDIF : 

Q9 = 0 

Q5 = AMINICTIME, 18.0) 


Wy 


NQAMNNMNNAANM AN 


YO 


75 


a rr gg mtg | 


DEXP(2*®Q5-6 .0) 
Q8/(Q8+.1/7Q8) 
Q7XCDEXPC(S1) 
RETURN 

END 


% 
a | 
ue te 


SUBROUTINE ALPHA 


THIS SETS UP THE INITIAL CONDITIONS OF THE 
PROBLEM INCLUDING THE DOMAIN AND THE 
ARTIFICIAL BOUNDARY 


THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER CIX = 20, IY = 80) 
PARAMETER CH = .025, DELTAT = .0325) 
PARAMETER CNDIM = 2*IX+1) 

THE FOLLOWING ARE COMPLEX ARRAYS 
COMPLEX*16 PC -IX: IX, O: IY, 3), WC -IX: IX, 3) 
COMPLEX LB, FACR 


COMPLEX AC -IX:1X, -IX: IX) 

COMPLEX AAC NDIM » NDIM ), BBC NDIM » NDIM 
COMPLEX ALC -40: 40) 

COMPLEX LAND 


COMMON/SS/P, W 
COMMON/TT/KP,KF,INFA, THETAI,L, X, PI, 
EISP, OMEGA, T 
COMMON/’CC/LB, FACR 
COMMON’ FF/A 
COMMON/GG/GAMMA, BETTA 
COMMON/HH/MINN, MAXX 
COMMON/II/AA, BB 
COMMON/KK/IPVT 
COMMON/RN/RUNS 
COMMON/FEL/U,V,TP 
COMMON’SN/FREQ, SPD 


REAL THETAI, INFA, KF, EISP, OMEGA, L, KP 

REAL XI, X, T 

REAL GAMMAC-40: 40), BETTA (€-40:40), PI’ 

REAL SPD, FREQ 

INTEGER I, U, V, TP, MINN, MAXX, IPVTCNDIM), RUNS, Y 
Ngee LY/2 

WRITEC2,75) 


FORMATC1X,' *) 
HEN = RUNSXDELTAT 


DO 10 M = MINN, MAXX,1 


ALCM) = CMPLX(0.0, 0.0) 
DO 30 I = -IX, IX, 1 


60 


) 


OM 


30 
10 
32 


111 


10 


1] 


30 


35 


G5 


XI = I*h 
IFC ABSCI) .EQ. IX) THEN 
DEL = .5 


ELSE 
DEL = 1.0 
ENDIF 
LAND = CEXP(-LBXYXHXBETTACM) + CLBXHEN)) 
AL(CM) = ALCM) + .5¥LANDXH*¥DELXCEXP(-LBXGAMMAC(M) ®XI)XPCI,Y,3) 
CONTINUE 
CONTINUE 
WRITEC2, 32) 
FORMATCIX, *ALPHA "N"S FOR RADIATING MODES‘) 
DO 111 M = MINN, MAXX 
WRITEC2,*®CABSCAL CM) ) 
CONTINUE 
RETURN 
END 


SUBROUTINE CLOSE 


THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 

PARAMETER CIX = 20, IY = 80) 
PARAMETER CH = .025, DELTAT = .0325) 
PARAMETER (CXYX= 40) 
PARAMETER CTNT= 1000) 
PARAMETER (CNDIM = 2X®IX+1) 

THE FOLLOWING ARE COMPLEX ARRAYS 
COMPLEXX16 PC -IX: IX, 0: ITY, 3), WO -IX: IX, 3) 
COMPLEX LB, FACR 

THE VALUES IN THE COMMON "BLOCKS" ARE VARIABLES 
COMMON/SS/P, W 


COMMON/TT/KP,KF,INFA, THETAI,L, X, PI, 
1 


EISP, OMEGA ,T 
COMMON/CC/LB, FACR 
COMMON/HH/MINN, MAXX 
COMMON/’GF/WGGRF , PGGRF 


REAL KP, KF, INFA, THETAI, L, PI, EISP, OMEGA 
REAL X, T 

REAL XX€ 0: G1), WPLOTC -IX: IX), 2Z(0:1000) 
REAL PPLOSC -IX: IX), SMALL(1000) 


INTEGER II, RUNS, SAM, IR 
PARAMETER (SAM = 1000) 


REAL WGGRF( SAM ), PGGRFC SAM ) 
DO 10 II = 0, XYX, 1 
XXCII) = II 
CONTINUE 
DO 11 IR = 0, TNT, 1 
ZZCIR) = IR 
CONTINUE 


DO 30 I = 1,1000,1 

SMALL(I) = WGGRFCI) 
CONTINUE 
DO 35 I = -IX, IX, 1 

PPLOSCI) = CDABS(P(I, IY,3)) 
CONTINUE 
DO 45 I = -IX, IX, 1 

WPLOTCI) = CDABS(WC(I,3)) 
CONTINUE 


6] 


AAaAAN: 


OOAMAAG O 


© 


OOOO 


34 


44 


88 


N= N = 


CALL COMPRS 

CALL PLOTD(ZZ,SMALL,1000, .TRUE. ,"LINLIN’,*CONV FACTORS’, 
‘CONVERGENCE GRAPHS$',*RUN NUMBER + 2000$', 
"CONVERGENCE VALUE$'*) 


CALL PLOTD(XX,PPLOS,41, .TRUE. ,‘LINLIN', FLUID PRESSURES', 


"PI72 LENGTH=1.2 FREQ= 3750$', HORIZONTAL POSITIONS", 
‘SCALED PRESSURES § ) 

CALL PLOTD(CXX,WPLOT,41, .TRUE. ,*LINLIN', *DISP AT 3750 HZ$', 
"SCALED LENGTH OF 1.2$',*SCALED HORIZONTAL POSITIONS', 
"SCALED DISPLACEMENT OF PLATES$') 

CALL DONEPL 

RETURN 

END 


SUBROUTINE FINN 


THIS SUBROUTINE MAKES SURE WE HAVE REACHED THE STEADY STATE 
SOLUTION OF THE PROBLEM, IE *EQUALIBRIUM' 
THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 


PARAMETER CIX = 20, IY = 80) 

PARAMETER CH = .025, DELTAT = .0325) 
THE FOLLOWING ARE COMPLEX ARRAYS 

COMPLEX*16 PC -IX: IX, 0: IY, 3), WC -IX: IX, 3) 

COMMON/SS/P, W 

COMMON/TT/KP,KF,INFA, THETAI, L, X, PI, 

EISP, OMEGA ,T 

COMMON/EE/WGTWNO, WGTHE, WGTOT, CON 

COMMON/PT/PGTWO, PGTHE, PGTOT 

COMMON/RN/RUNS . 

COMMON/GF/WGGRF , PGGRF 


eae KP, L, OMEGA, EISP, KF, INFA, PI, THETAI 
EA X» 

REALX8 WGTWO,WGTHE,+HGTOT,CON 

REALX&8 PGTWO,PGTHE,?PGTOT 

INTEGER I, RUNS, SAM 

PARAMETER CSAM = 1000) 

REAL WGGRFC SAM 3, PGGRFC SAM ) 


DO 34 1 = -IX+2, IX-2, 1 
WGTHO = WCI,2)*DCONJGCHCI,2))+ WGTHO 
WGTHE = WCI,3)XDCONJGCWCI,3))+ WGTHE 
WGTOT =  (WGTHE - WGTWO) 7 WGTHE 

CONTINUE 

DOe 44 1 => ae Doe 
PGTHO = PCI,IY,2)*DCONJG(PCI,1Y,2))+ PGTWO 
PGTHE = PCI,IY,3)*XDCONJGCP(I,IY,3))+ PGTHE 
PGTOT = | (PGTHE - PGTWO) 7 PGTHE 

CONTINUE 


IF CRUNS.GT. 1999) THEN 


WGGRF CRUNS-1999) = WGTOT 
PGGRFCRUNS-1999) = PGTOT 
ENDIF 

WRITEC2,88) 


FORMATC1X, *‘WGGRF,PGGRF IN FINN®) 
WRITEC2,%) WGGRFCRUNS-1999) 
WRITEC2,%) PGGRFCRUNS-1999) 
RETURN 
END 
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